/*
 * MandelbrotQuadTreeCompute.c
 *
 *  Created on: 15.07.2012
 *      Author: 116
 */

#include<stdlib.h>
#include <stdio.h>

#include "AreaComputationInformation.h"
#include "MandelbrotQuadTree.h"
#include "MandelbrotQuadTreeHelpers.h"

MandelbrotQuadTree * nextNodeToCompute(MandelbrotQuadTree * sourceNode, unsigned int maxComputationDepth)
{
    //id the source node hasn't been computed : compute it
    if (sourceNode->status == VOID)
    {
        return sourceNode;
    }
    // we want to continue computation only on nodes which are not completely outside or inside the mandelbrot set
    else if (sourceNode->status == BROWSED)
    {
        //if the node has no children : create them
        if (!sourceNode->children)
        {
            splitNode(sourceNode);
            MandelbrotQuadTree * nextNode = &(sourceNode->children[TopLeft]);
            if (getNodeDepth(nextNode) > maxComputationDepth)
            {
                return 0;
            }
            else
            {
                return nextNode;
            }
        }
        else
        {
            MandelbrotQuadTree * n1 = nextNodeToCompute(&(sourceNode->children[TopLeft]), maxComputationDepth);
            MandelbrotQuadTree * n2 = nextNodeToCompute(&(sourceNode->children[TopRight]), maxComputationDepth);
            MandelbrotQuadTree * n3 = nextNodeToCompute(&(sourceNode->children[BottomLeft]), maxComputationDepth);
            MandelbrotQuadTree * n4 = nextNodeToCompute(&(sourceNode->children[BottomRight]), maxComputationDepth);

            MandelbrotQuadTree * best = 0;

            //try to assign at least 1 non NULL pointer
            if (n1)
            {
                best = n1;
            }
            else if (n2)
            {
                best = n2;
            }
            else if (n3)
            {
                best = n3;
            }
            else if (n4)
            {
                best = n4;
            }

            if (!best)
            {
                // if we didn't get any non NULL pointer, then all hope is lost :( --> return NULL
                return 0;
            }

            //get the least deep node : breadth first browsing
            //don't test for n1 because is best!=n1, then n1 is NULL
            if (n2 && getNodeDepth(n2) < getNodeDepth(best))
            {
                best = n2;
            }
            if (n3 && getNodeDepth(n3) < getNodeDepth(best))
            {
                best = n3;
            }
            if (n4 && getNodeDepth(n4) < getNodeDepth(best))
            {
                best = n4;
            }

            if (getNodeDepth(best) > maxComputationDepth)
            {
                return 0;
            }
            else
            {
                return best;
            }
        }

        return 0;
    }
    else
    {
        return 0;
    }
}

void borderPointsAsDouble(MandelbrotQuadTree * node, double * destArray, unsigned int pointsPerSide)
{
    double minX = node->minX;
    double maxX = node->maxX;
    double minY = node->minY;
    double maxY = node->maxY;

    double step = (maxX - minX) / (double) (pointsPerSide - 1);

    // bottom side of the rectangle
    unsigned int baseIndex = 0;
    for (unsigned int i = 0; i < pointsPerSide; i++)
    {
        destArray[baseIndex + 2 * i] = (minX + (double) i * step);
        destArray[baseIndex + 2 * i + 1] = (minY);
    }
    baseIndex += 2 * pointsPerSide;

    // top side of the rectangle
    for (unsigned int i = 0; i < pointsPerSide; ++i)
    {
        destArray[baseIndex + 2 * i] = (minX + (double) i * step);
        destArray[baseIndex + 2 * i + 1] = (maxY);
    }
    baseIndex += 2 * pointsPerSide;

    // left side of the rectangle
    for (unsigned int i = 0; i < pointsPerSide; ++i)
    {
        destArray[baseIndex + 2 * i] = (minX);
        destArray[baseIndex + 2 * i + 1] = (minY + (double) i * step);
    }
    baseIndex += 2 * pointsPerSide;

    // bottom side of the rectangle
    for (unsigned int i = 0; i < pointsPerSide; ++i)
    {
        destArray[baseIndex + 2 * i] = (maxX);
        destArray[baseIndex + 2 * i + 1] = (minY + (double) i * step);
    }
}

int isInsideMandelbrotSet(double * array, unsigned int pointsCount, unsigned int maxIter)
{
    for (unsigned int i = 0; i < pointsCount; ++i)
    {
        double real = array[2 * i];
        double img = array[2 * i + 1];
        double realsqr = real * real;
        double imgsqr = img * img;

        double real1 = real;
        double img1 = img;
        double real2, img2;

        unsigned int iter = 0;
        while ((iter < maxIter) && ((realsqr + imgsqr) < 4))
        {
            real2 = real1 * real1 - img1 * img1 + real;
            img2 = 2 * real1 * img1 + img;

            real1 = real2 * real2 - img2 * img2 + real;
            img1 = 2 * real2 * img2 + img;

            realsqr = real2 * real2;
            imgsqr = img2 * img2;
            real1 = realsqr - imgsqr + real;
            img1 = 2 * real2 * img2 + img;

            iter += 2;
        }

        if (iter < maxIter)
        {
            return 0;
        }
    }

    return 1;
}

void innerPointsAsDouble(MandelbrotQuadTree * node, double * array, unsigned int pointsPerSide)
{
    double minX = node->minX;
    double maxX = node->maxX;
    double minY = node->minY;
    double maxY = node->maxY;

    double stepX = (maxX - minX) / (double) (pointsPerSide - 1);
    double stepY = (maxY - minY) / (double) (pointsPerSide - 1);

    // bottom side of the rectangle
    for (unsigned int i = 0; i < pointsPerSide; i++)
    {
        int base = i * pointsPerSide;
        double xVal = minX + (double) i * stepX;
        for (unsigned int j = 0; j < pointsPerSide; ++j)
        {
            array[2 * (base + j)] = xVal;
            array[2 * (base + j) + 1] = (minY + (double) j * stepY);
        }
    }
}

AreaComputationInformation isOutsideMandelbrotSet(double * array, unsigned int pointsCount, unsigned int maxIter, unsigned int differentIterLimit)
{
    unsigned int min, max;
    //init min and max iter
    {
        double real = array[0];
        double img = array[1];
        double realsqr = real * real;
        double imgsqr = img * img;

        double real1 = real;
        double img1 = img;
        double real2, img2;

        unsigned int iter = 0;
        while ((iter < maxIter) && ((realsqr + imgsqr) < 4))
        {
            real2 = real1 * real1 - img1 * img1 + real;
            img2 = 2 * real1 * img1 + img;

            real1 = real2 * real2 - img2 * img2 + real;
            img1 = 2 * real2 * img2 + img;

            realsqr = real2 * real2;
            imgsqr = img2 * img2;
            real1 = realsqr - imgsqr + real;
            img1 = 2 * real2 * img2 + img;

            iter += 2;
        }
        min = iter;
        max = iter;
    }

    for (unsigned int i = 0; i < pointsCount; ++i)
    {
        double real = array[2 * i];
        double img = array[2 * i + 1];
        double realsqr = real * real;
        double imgsqr = img * img;

        double real1 = real;
        double img1 = img;
        double real2, img2;

        unsigned int iter = 0;
        while ((iter < maxIter) && ((realsqr + imgsqr) < 4))
        {
            real2 = real1 * real1 - img1 * img1 + real;
            img2 = 2 * real1 * img1 + img;

            real1 = real2 * real2 - img2 * img2 + real;
            img1 = 2 * real2 * img2 + img;

            realsqr = real2 * real2;
            imgsqr = img2 * img2;
            real1 = realsqr - imgsqr + real;
            img1 = 2 * real2 * img2 + img;

            iter += 2;
        }

        if (iter < min)
        {
            min = iter;
        }
        else if (iter > max)
        {
            max = iter;
        }

        if ((max - min + 1) > differentIterLimit)
        {
            AreaComputationInformation area;
            area.maxIter = -1;
            area.minIter = -1;
            return area;
        }
    }

    AreaComputationInformation area;
    area.maxIter = max;
    area.minIter = min;

    return area;
}

int computeNextNode(MandelbrotQuadTree * root, double * edgeArray, double * innerArray, unsigned int pointsPerSide, unsigned int maxComputationDepth,
        unsigned int maxIter, unsigned int diffIterLimit)
{
    MandelbrotQuadTree * nextNode = nextNodeToCompute(root, maxComputationDepth);
    if (nextNode)
    {
//        printf("\nNext node:\n");
//        show(nextNode);

        borderPointsAsDouble(nextNode, edgeArray, pointsPerSide);

//        printf("Got edge points\n");

        int inside = isInsideMandelbrotSet(edgeArray, pointsPerSide * 4, maxIter);

//        printf("Inside=%i\n", inside);

        if (!inside)
        {
            innerPointsAsDouble(nextNode, innerArray, pointsPerSide);

//            printf("Got inner points\n");

            AreaComputationInformation area = isOutsideMandelbrotSet(innerArray, pointsPerSide * pointsPerSide, maxIter, diffIterLimit);

            if (area.maxIter != -1)
            {
//                printf("Outside: min=%i, max=%i\n", area.minIter, area.maxIter);
                nextNode->status = OUTSIDE;
                nextNode->minIter = area.minIter;
                nextNode->maxIter = area.maxIter;
            }
            else
            {
//                printf("Browsed\n");
                nextNode->status = BROWSED;
            }
        }
        else
        {
            nextNode->status = INSIDE;
        }
        return getNodeDepth(nextNode);
    }
    else
    {
        printf("No next node\n");
        return -1;
    }
}
